** Spatial approach, 3: Define Neighborhood ** 

		**********------------------------------**********
	* Contents: 
		* Build weight matrices for a given neighborhood definition
		* Grid corresponds to "Step function" in final paper  
		* Input: cell-level data based on spatial_2_collapse_s 
		* Output: weight matrixes in .dta format 

capture set maxvar 30000, permanently // required to run file.  
	
capture log close 
clear all
set more off 

cd "R:\WSV2\TBu_AKe\Spatial_NEW"

capture mkdir Data
global store "R:\WSV2\TBu_AKe\Spatial_NEW\Data"

log using spatial_3_wmat_grid, replace 


cd "R:\WSV2\TBu_AKe\Spatial_NEW"


use "C:\Users\hy65byfe\Desktop\smerge_0712\spatial_2_collapse_s.dta" // data set with 4950 cells, created in spatial_2_*  


scalar p = 0 // if 0, straight distance without decay : x^0 = 1, so all are equal.
scalar sigma = 47*0.59 // sigma at comp. line 
scalar kappa = 51.7*0.59 // kappa at comp. line 
scalar D = 3 // width of increment 

***************************
**** reduce dimensions ****
***************************

keep if year == 2011 & ccode == 1 // reduce data in memory 
  
sort cell // now sorted: first bunching (1-50), then restricted (51-1015) then rest. 

gsort -bunching +cell 
gen bcell = _n if bunching == 1 // identify B

gsort -restricted +cell 
gen rcell = _n if restricted == 1 // identify R


**********************************************
*** define upper and lower bounds of cells ***
**********************************************
		
	** gen mid-point for E 
*gen a_mid = a_bin 
*gen e_mid = `=kappa' + `=sigma'*a_mid //  

	** gen upper and lower level of E at mid-point (e.g. 6kg, 6.5kg ...)
replace e_bin = k_bin + a_bin*`=sigma' // gravity center: upper edge of mid-point 
gen e_ml = k_bin + a_bin*`=sigma' - `=D' // lower edge of mid-point 


replace eei_bin = (e_bin/(51.7+47*a_bin))* 100 // assert that bins are consistent with eei boundaries. 

	** label as safety net 
label variable e_bin "E_i at midpoint"

label variable a_min "a_i left"
label variable a_max "a_i right"
label variable k_bin "kappa (std. E)"


keep eei_bin a_* e_* k_bin s_bin eei_bin cell bcell rcell // cut 

sort cell // order is crucial 
xtset cell 
spset cell 

***********************************
******* start loops ***************
***********************************

************************************	 
*  3a: neighborhood distance count *
************************************

gen k_run = (s_bin/3) // running index with 0,1,2,3,4,5,6 ... for cells in E-direction 

gen a_run = (a_bin*2 - 6) +1 // running index with 1,2,3, ... for cells in a-direction (must start at 1)

* scalar list 

**************************************
* 3b: define scalar for every cell ***
**************************************


sum cell, meanonly 
scalar cmax=`r(max)'
local c = `r(max)' // scalar for each cell.  
	
sum bcell, meanonly 
scalar bmax=`r(max)'

sum rcell, meanonly 
scalar rmax=`r(max)'

cd "$store"

forvalues q = 5(5)40 {
	preserve 

	scalar q = `q' // after q cells, jump to next strip 

	sum k_run, meanonly
	scalar max_k = `r(max)'+`=q' // 
	scalar min_k = `r(min)'

	sum a_run, meanonly
	scalar max_a = `r(max)' +`=q' // going left after 1to1 in grid 
	scalar min_a = `r(min)'

	*** create group address *** 

	egen k_run2 =cut(k_run),at(0(`=q')`=max_k') // first 5 addresses get 1, second 5 get 2: like city blocks. 

	egen k_group = group(k_run2) 
	replace k_group = k_group - 1 // start at zero 

	*sum k_run if missing(k_group) // all negative. 


**************************************
* 3b: define scalars for every cell ***
*************************************

	forvalues i = 1/`=cmax' {
 	
		* a-address  
		quietly sum a_run if cell == `i'
		local a = r(mean)
		scalar a`i' = `a'
	
		* std.  
		quietly sum k_group if cell == `i'
		local k = r(mean)
		scalar k`i' = `k'
	
		* eei 
		quietly sum eei_bin if cell == `i'
		local x = r(mean)
		scalar x`i' = `x'

	}

	
	*************************
	* 4a: define nieghbors ***
	*************************

	forvalues m = 1/`=bmax' {
     quietly sum cell if bcell==`m' //  
     local i = r(mean)
	  * i contains cellid of cell i in B 	 
          forvalues l = 1/`=cmax' {
          quietly sum cell if cell==`l'
          local j = r(mean)
		  * j contains cellid of cell j in R 
		  
		  * knock-out condition 
		  scalar x_`i'_`j' 	= cond( x`j' > x`i'   &  inrange(x`j', 59, 87)		,	1,0) // only neighbors in B. 
		 
		  * define neighbors 
		  scalar no_`i'_`j' =  cond( -1*(a`i' - a`j') 	== k`j',       1,0) // nord-ost : i < j 
		  scalar nw_`i'_`j' =  cond(    (a`i' - a`j') 	== k`j',       1,0) // nord-west: i > j  
		  
		  }
        
     }	 


	*********************************
	*** 4b: transform to matrix *****
	*********************************

	*** neighborhood matrices: MANUAL: Sektors 
	mat define X  = I(`=cmax') // knock out 
	mat define O`q' = I(`=cmax') // east  
	mat define W`q' = I(`=cmax') // sw 2


	*********************
	**  5: fill matrix **
	*********************


	forvalues m = 1/`=bmax' {
		quietly sum cell if bcell==`m' & bcell != . 
		local i = r(mean)
		forvalues l = 1/`=cmax' {
			quietly sum cell if cell==`l' 
			local j = r(mean)
		  ** run through combinations and put in distance **
			mat X[`m',`l']=(x_`i'_`j')
			mat O`q'[`m',`l']=(no_`i'_`j')
			mat W`q'[`m',`l']=(nw_`i'_`j')
         }
		forvalues c = 1/`=cmax' {
		** replace own-cross with zero (identical cell)
		mat X[`c',`c']=0
		mat O`q'[`c',`c']=0
		mat W`q'[`c',`c']=0
			}
    } 
	
	matrix XO`q'=hadamard(X, O`q')
	matrix XW`q'=hadamard(X, W`q')


	svmat XO`q', names(o`q'_)
	svmat XW`q', names(w`q'_)

	spmatrix fromdata O`q' = o`q'_1 - o`q'_4950, normalize(none)
	spmatrix save O`q' using O`q'.stswm, replace

	spmatrix fromdata W`q' = w`q'_1 - w`q'_4950, normalize(none)
	spmatrix save W`q' using W`q'.stswm, replace

	save spatial_3_wmat_grid_`q', replace

	drop k_run2
	drop k_group 
	
	restore 

}


********************** repeat for horizontal and vertical matrix *******************************

	
	scalar q = 0 // no jumps  
	local q = 0 //  local macro 

	sum k_run, meanonly
	scalar max_k = `r(max)'+`=q' // 
	scalar min_k = `r(min)'

	sum a_run, meanonly
	scalar max_a = `r(max)' +`=q' // irrelevant for this case, q=0. 
	scalar min_a = `r(min)'

	
	**************************************
	* 3b: define scalars for every cell ***
	*************************************


	forvalues i = 1/`=cmax' {
 	
		* a-address  
		quietly sum a_run if cell == `i'
		local a = r(mean)
		scalar a`i' = `a'
		
		/* 
		* std.  
		quietly sum k_group if cell == `i'
		local k = r(mean)
		scalar k`i' = `k'
		*/ 
		
		* eei 
		quietly sum eei_bin if cell == `i'
		local x = r(mean)
		scalar x`i' = `x'
		
		* e_bin (low-left corner)
		quietly sum e_bin if cell == `i'
		local e = r(mean)
		scalar e`i' = `e'

	}


	forvalues m = 1/`=bmax' {
     quietly sum cell if bcell==`m' // 
     local i = r(mean)
	  * i contains cellid of cell i in B 	 
          forvalues l = 1/`=cmax' {
          quietly sum cell if cell==`l'
          local j = r(mean)
		  * j contains cellid of cell j in R 
		  
		  * knock-out condition  
		  scalar x_`i'_`j' 	= cond( x`j' > x`i'   &  inrange(x`j', 59, 87)		,	1,0) // only neighbors in B. 
		 
		  
		  * define neighbors 
		  scalar r_`i'_`j' =  cond(a`i' ==  a`j',			       1,0) // rook 
		  scalar h_`i'_`j' =  cond(abs(e`i'-e`j') <= 6,       1,0) //    
		  }
        
     }	 


	*********************************
	*** 4b: transform to matrix *****
	*********************************

	*** neighborhood matrices: 
	mat define X  = I(`=cmax') // knock out 
	mat define R = I(`=cmax') //  vertical / rook 
	mat define H = I(`=cmax') // horizontal


	************************
	**  5: fill matrix *****
	************************

	forvalues m = 1/`=bmax' {
		quietly sum cell if bcell==`m' & bcell != . 
		local i = r(mean)
		forvalues l = 1/`=cmax' {
			quietly sum cell if cell==`l' 
			local j = r(mean)
		  ** run through combinations and put in distance **
			mat X[`m',`l']=(x_`i'_`j')
			mat R[`m',`l']=(r_`i'_`j')
			mat H[`m',`l']=(h_`i'_`j')
         }
		forvalues c = 1/`=cmax' {
		** replace own-cross with zero (identical cell)
		mat X[`c',`c']=0
		mat R[`c',`c']=0
		mat H[`c',`c']=0
			}
    } 
	
	matrix XR=hadamard(X, R)
	matrix XH=hadamard(X, H)


	svmat XR, names(r_)
	svmat XH, names(h_)

	spmatrix fromdata R = r_1 - r_4950, normalize(none)
	spmatrix save R using R.stswm, replace

	spmatrix fromdata H = h_1 - h_4950, normalize(none)
	spmatrix save H using H.stswm, replace

	save spatial_3_wmat_grid_polar, replace

	
log close 

clear 

exit 